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A Cartesian grid method combined with a simplified gas kinetic scheme is presented for subsonic 
and supersonic viscous flow simulation on complex geometries. Under the Cartesian mesh, the 
computational grid points are classified into four different categories, the fluid point, the solid point, 
the drop point, and the interpolation point. The boundaries are represented by a set of direction- 
oriented boundary points. A constrained weighted least square method is employed to evaluate the 
physical quantities at the interpolation points. Different boundary conditions, including isothermal 
boundary, adiabatic boundary, and Euler slip boundary, are presented by different interpolation 
strategies. We also propose a simplified gas kinetic scheme as the flux solver for both subsonic and 
supersonic flow computations. The methodology of constructing a simplified kinetic flux function 
can be extended to other flow systems. A few numerical examples are used to validate the Cartesian 
grid method and the simplified flux function. The reconstruction scheme for recovering the boundary 
conditions of compressible viscous and heat conducting flow with a Cartesian mesh can provide a 
smooth distribution of physical quantities at solid boundary, and present an accurate solution for 
the flow study with complex geometry. 


I. INTRODUCTION 

Solving flow problem on conformal mesh or unstructured mesh is the most practical strategy in industry 
applications. However, the grid generation on complex geometry is time consuming and requires sophisti¬ 
cated technique. In this sense, the grid generation on complex computational domain is still the bottle neck 
of practical applications. Thereby, the Cartesian grid method arises, in which the flow region is discretized 
by a Cartesian grid regardless of the shapes of objects inside the flow region. The obvious advantage of 
this method over the conventional conformal approach is that the computational mesh is easily generated 
when different geometries are considered. Cartesian grid methods free the researchers and engineers from 
the burdensome grid generation, but introduce two new problems about boundary treatment. 

The first problem with the Cartesian grid is about how to represent the boundary. It depends on the 
boundary property. For the multi-fluid and multi-phase flows, where the boundaries or interfaces are de¬ 
formable, the volume-of-fluid (VOF) method^ and the phase-field approach^ are the most popular ap¬ 
proaches. The boundary is reconstructed from an auxiliary function, say, a marker function. Or, the 
boundary is captured by solving a partial differential equation for the phase field. These approaches are 
efficient for boundary deformation and splitting. But obviously, this kind of approaches is less accurate than 
the Lagrangian representation of the boundary. The latter approach takes the boundary as a sharp interface 
either explicitly tracking as curves^ or as level sets^. Sharp interface is desirable for high Reynolds number 
flows where the boundary layer plays an important role. 

The second problem is how to impose the boundary condition. Following the classification by Mittal and 
laccarino^, there are two different categories. The first one is the continuous forcing approach; the second 
one is the discrete forcing approach. The continuous forcing approach directly modifies the governing 
equation by adding a forcing term to take the boundary effect into account^ii. An ideal force term is 
represented by a Dirac delta function. Since the boundary cuts the grid line at arbitrary location, the 
forcing should be distributed over a band of cells around the boundary point. This approach results in a 
diffusive boundary. The discrete forcing approach imposes the boundary condition on the numerical solution 
directly at discrete level. Mittal and laccarino^ further subdivided the second category into ’’Indirect BC 
Imposition” and ’’Direct BC Imposition”. The former one employs a forcing term which is determined from 
a priori estimation of flow field at discrete level^. The external force is explicitly computed in advance 
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or solved by an implicit method to guarantee the no-slip boundary conditioned^. However, the artificial 
forcing procedure diffuses the flow field around the boundary. For this reason, the ’’Direct BC Imposition” 
approaches retain the boundary as a sharp interface with no spreading. This can usually be accomplished 
by modifying the computational stencil near the boundary to directly impose the boundary condition. This 
kind of approach is always referred to as Cartesian grid method which is the main focus in this study. 

Berger and LeVequeii presented a Cartesian mesh algorithm with adaptive refinement to compute flows 
around arbitrary geometries by solving the Euler equations. They treated the intersection between the 
grid line and the boundaries as a grid point and performed conventional finite difference scheme with 
amendments of boundary fluxes. However, the small cell instability was observed. Pember et al^ proposed a 
corrector applied to cells at the boundary to redistribute flow quantities in order to maintain the conservation, 
therefore, avoid time step restrictions arising from small cells. In these methods, the amendments or the 
correctors depend on the specific physical system, thereby, the numerical scheme is difficult for further 
extension. 

Forrer and Jeltscbi^ proposed a ghost cell method to apply the boundary condition on Cartesian grid for 
inviscid compressible Euler equations. The flow quantities at ghost cells which locate inside the boundary 
are set by their images in the fluid side across the boundary. As a result, the grid points near the boundary 
also preserve a complete control volume. Thus the small cell restriction is removed. Hereafter, this kind of 
method is referred to as the boundary interpolation method. In recent years, the boundary interpolation 
method receives more and more attentions. Udaykumar et al^ proposed a finite-difference formulation 
to track solid-liquid boundaries on a fixed underlying grid. The interface is not of finite thickness but is 
treated as discontinuity and is explicitly tracked. The imposition of boundary conditions exactly on a sharp 
interface is performed using simple stencil readjustments in the vicinity of the interface. Ghias et alr^ 
proposed an immersed boundary method for computing viscous, subsonic compressible flows with complex 
shaped stationary immersed boundaries. The method also employs a ghost cell technique for imposing the 
boundary conditions on the immersed boundaries. As Gibou et aL— pointed out that the single directional 
interpolation scheme is poorly behaved for small cells. Therefore, Ghias et al. turned to a multidimensional 
bilinear interpolation and evaluated the flow variables at ghost cell to avoid the small cell problem. 

Instead of interpolating variables to the ghost cell which is located inside the solid, some researchers 
interpolate flow variables at the fluid points, or, the interface points which are very close to the boundary 
at the fluid side. Yei^ et al. and Udaykumar et alM- proposed a sharp interface Gartesian grid method 
for solving the incompressible Navier-Stokes equations with complex moving boundaries. The boundaries 
cut the mesh, then form some trapezoidal cells. If the trapezoidal cell is very small, then it merges with 
its neighboring cell, and forms an irregular control volume. The finite volume method can be applied on 
the merging control volume. Then the small cell problem is circumvented. They proposed an interpolation 
procedure to construct the flow field near the boundary, so that the boundary condition is satisfied. Tullio 
et alr^ and Palma et alr^ combines the method for solving the three-dimensional preconditioned Navier- 
Stokes equations for compressible flows with an immersed boundary approach, to provide a Gartesian-grid 
method for computing complex flows over a wide range of the Mach number. Moreover, a flexible local 
grid refinement technique is employed to achieve high resolution near the immersed body and in other 
high-flow-gradient regions. 

Actually, many Gartesian grid method solving the compressible Navier-Stokes equations are based on the 
interpolation method to handle the boundary conditions. Peng et al^ has compared the forcing method 
and the interpolation method in the framework of lattice Boltzmann method. They concluded that the 
interpolation method perform better than the forcing method. In recent years, more delicate interpolation 
procedures are employed for boundary treatment. Seo and Mittal^i proposed a new sharp-interface immersed 
boundary method based approach for the computation of low-Mach number flow-induced sound around 
complex geometries. Their method applies the boundary condition on the immersed boundary to a high 
order by combining the ghost cell approach with a weighted least square error method based on a high order 
approximating polynomial. Duan et al^ extended the Cartesian grid method to high order accuracy, and 
simulated turbulent flow at high Mach numbers. 

In this paper, we are going to further develop the boundary interpolation method on Cartesian grid, and 
combines it with well developed gas kinetic schemes. The main purposes are presented as follows: 

1. simplify the gas kinetic schemes, and reduces the computational cost. 

2. formulate a constrained weighted least square^^ interpolation for different boundary condition. 
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II. SIMPLIFIED GAS KINETIC SCHEME 


Gas kinetic schemes (GKS) were developed by Xu et It is based on kinetic theory, using particle 

distribution function to simulate macroscopic flow motion. Gas kinetic scheme can be considered as a flux 
solver of a generalized Riemann problem. In this study, the Bhatnagar-Gross-Krook kinetic equation^ is 
employed. It takes the following form. 


dt 


^ 

5x r ’ 


( 1 ) 


where /(x, t, u, represents the particle velocity distribution function. It is a function of location x, particle 
velocity u, the internal degree of freedom ^ and time t. The right hand side of Eq.(IT]) represents the relaxation 
process, g is the equilibrium state which reads, 


K + 2 


( 2 ) 


where p denotes density, U represents the macroscopic velocity, T denotes temperature, R is gas constant 
and K is the equivalent internal degree of freedom which is 1 for monatomic gas, and 3 for diatomic gas in 
two dimensional simulation. The conservative variables W can be derived by taking moments of /, 


W = 



= J -tpldud^ =< f > 


where tp = (l,u, and E denotes the energy per unit mass. The macroscopic flux of a given 

distribution function / can be written as. 


F = 


'0n • ufdud^ =< n • u/ >, 


where n denotes a certain direction, say, the normal direction of a surface. 

The original GKS is a finite volume method. If the numerical flux F is known at the cell interfaces of a 
cell, the macroscopic variables for the next time step becomes. 


= W” 



( 3 ) 


where Sk denotes the area of the fcth surface of the cell, and V is the volume of the cell. In this study, we 
adopt the finite difference formulation. Then, for two dimensional problems, 

= wr, - (4) 

^ 2 ; Ay 

Figure [T] shows the stencil for the numerical flux Fj_|_]^/ 2 j- The solid square represents the flux point, or the 
cell interface in a finite volume method. To calculate the numerical flux, the distribution function at the flux 
point is needed. The original GKS employs the local analytical solution of the Eq. o as the distribution 
function, which is shown as follows, 

/(xo, t, u, = e"*/^/(xo - ut, t, u, 

+ - [ (5) 

T Jo 

where x' = Xq — u{t — t'). This solution includes the complete information from both space and time. 
The numerical method based on this solution performs well in both continuous and discontinuous flows^^. 
However, the numerical cost is a little higher in comparison with classical Riemann solver-based finite volume 
method. And also, a simplified GKS was proposed for low speed and continuous flow simulation in which 
the Ghapman-Enskog expansion is directly employed as the local solution to evaluate the numerical flux at 
the cell interface^^. 


fit) = 9c- r(u • gx + gt) + 5 x • X -k gtt. 


(6) 
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FIG. 1. The stencil for numerical flux in x direction. 


The variables Xq, u, ^ are omitted hereafter for the sake of simplicity, gc is the equilibrium state at the cell 
interface. It is constructed by central difference scheme, namely gc = A4[Wc]. The subscript c denotes the 
central difference scheme. This continuous reconstruction is accurate for low speed flow simulation where 
no shock exists. But, it cannot be used in supersonic flow simulation. 

In this study, we replace the leading order term gc by a non-equilibrium distribution function /o, 

f{t) = fo- t(u • gx + 5 t) + 5 x • X -t- gtt. (7) 

The first term on the right hand side represents the upwind effect, the second term represents the viscous 
effect, the third term represents the time evolution up to second order accuracy. Similar scheme has been 
proposed^i, in which the leading order term is (1 — a)fo + age- Figure [T] shows the stencil used in 2D 
reconstruction. A uniform mesh is assumed. The non-equilibrium state /o is 


(gi=M[Wi], u>0 
= AI[Wr], n < 0 ’ 


( 8 ) 


where the subscript T and ’r’ represent the left and right side of the flux point respectively. W; and 
are constructed along the grid line. For the fluid side, the 3rd order WENO reconstruction is adopted, 


Wi = 


w-iW-i -\- wobFo 

W-i Wq 
WqWq + WiWi 


Wr = 

Wq -I- Wi 

where W represents every single macroscopic quantity of W. The weights are given as follows, 

1 

w_i = 


iisU + sY 


^0 A/ 2 \ ’ 

4(sf e) 


Wi = 


4(sf+i+e)’ 


where e is a small positive number to avoid zero denominator, and Si = Wi+ij — Wij, 

W-i = 

kFo = 2^i,j + 2 ^*+10’ 

Wi = - ^W.+2,,. 


( 9 ) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 

(16) 
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For the stencil contains drop point or interpolation point near the boundary, the van-Leer slope limiter is 
adopted to construct the flow quantities at either side of cell interface. The spatial derivative in the second 
term of Eq. 0 can be expressed in the expansion of equilibrium state, which formally reads, 

5x = a.^Af[W"], (17) 


where the expansion coefficient a reads. 


a-ip = ao,d + ai,d • u + a 2 ,d 



d = l,...,D. 


(18) 


where ai^^ = ■■■■, cii^du), and D denotes the dimension. For example, for two dimensional problems, 

D = 2. Then take the conservative moments of the spatial derivatives of equilibrium state, we have. 


a- < V'MiW”] >= W,,, 


(19) 


where W" denotes the conservative variables at the beginning of the time step. So, the corresponding 
conservative variables are W” =< fo > for the Eq. O- Then Wx is derived by a central difference scheme. 
Take two dimensional problems as an example. 




w —_ 

^ " 4Ay 


Ax 

(W.., 


i+i 


W 


i+i.i+i 


-W„_i -W,+i,,_i) 


( 20 ) 

( 21 ) 


Solving the Eq. (HH), we get the expression of g^- 

The temporal derivative is estimated by the direct difference of the equilibrium state on the time interval. 




( 22 ) 


where denotes the conservative variables at the end of the time step. It can be derived by applying 

the conservation constraint, < T{gt + u • g^) >= 0 at a cell interface, that is, 

W"+i = W” - At < u ■ 5x > ■ (23) 


The time dependent distribution function at flux point is, 

-7W[W"], Ad[W"+i] - Ad[W" 
fi.t) = fo-T{u-g^ + -1-- ‘-) + t - 


At 


At 


Therefore, the numerical flux reads. 


Fdt = n • 


ipuf {t)dud^dt 


= n • {At < u/o > —rAt < uu • t/x > 

+ (^ - ^) < u(M[W"+i] - M[W"]) >), 


(24) 


Up to this point, the simplified gas kinetic scheme is completed. Similar scheme has been studied in^^. 


III. CARTESIAN GRID METHOD 

Figure [2] illustrates the stencil for the boundary interpolation. The circles represent the fluid points; the 
squares represent the drop points which are excluded from the flux stencil and the interpolation stencil; the 
diamonds represent the solid points, while the empty diamonds are the interpolation points or ghost points; 
and the boundary points are denoted by points with their normal directions. The right top part of figure 
[2] illustrates how to classify a grid point. The grid points (0,D,C) are solid points, since they locate at 
the back of the neighboring boundary points (E and F). The grid points A and B are fluid points, which 
locate at the front of boundary points (E or F). If any grid point is too close to the boundary point, then, 
this point is marked as a drop point. When a solid point is involved in a flux stencil, then it becomes 
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FIG. 2. The stencil for the boundary interpolation. 


an interpolation point. After defining different kind of grid points, all the flux stencils are composed of 
interpolation points and fluid points. Meanwhile, the small cells formed by the cutting process are totally 
avoided, and no stability problem is observed. 

The boundary condition is satisfied by assigning the corresponding values at the interpolation points or 
ghost points. The solid points which are involved in a flux stencil are employed as interpolation points. It 
means that an interpolation point can preserve different values when it is in different flux stencil. To satisfy 
the boundary condition, the fluid variables at the interpolation point should be assigned by appropriate 
interpolation. As shown in figure [21 we use the fluid points inside the big circle to reconstruct the fluid 
variables at interpolation points to enforce the corresponding boundary conditions. 


A. Constrained weighted least square interpolation 


Seo^i proposed a weighted least square interpolation for low Mach number acoustic simulation. For low 
Mach number flow, the boundary value can be predetermined, since the density is constant, the velocity is 
given at boundary, and the pressure is totally determined by a Poisson equation. This is not the case for 
the compressible viscous flow simulation. The boundary variables are unknown before the interpolation. As 
a result, the weighted least square method cannot be used in compressible flow simulation directly. The 
boundary variables should be predetermined from the neighboring fluid properties. For example, Palma 
et assumed the constant pressure across the first layer of fluid point near the boundary, then use the 
derived boundary variables to construct the fluid variables at the interpolation point. In this study, we 
employ a constrained weighted least square method^^ to construct the fluid variables at the interpolation 
point. This interpolation process interpolates an intermediate flow state first, which is subsequently modified 
in order to satisfy the boundary condition. Then, the distribution of flow variables around the boundary is 
determined. Assume that the fluid variable around the boundary can be expressed as a polynomial, 

P{x,y)= X! > 0,9 > 0, (25) 

p+g<m 

where m is the order of polynomial. Given the values {ui} are known at a set of locations {{xi,yi)\i = 
I,N, N > dim(X)}. And they should be approximately satisfied in the sense of minimizing a weighted 
mean square error function, which reads. 


N 


minimize 




2=1 


where Wi is the weight which is written in an exponential form, 

(3'^iixi - xb)'^ + (jji - ytY) 


), 


(26) 


Wi = exp( 


(max(Aa;, A?/))^ 


(27) 
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where f3 is an adjustable parameter. And Vi denotes the every single variable of W^. 


X = (co,o, 

A, = (l, .. 


5 ^0,li 

^iViy 


Cm—1 

r^rn-l 


Vi-. 



At first, the unconstrained problem described above is solved. Let v = (vi,'y 2 , 


(28) 

(29) 



/ Wo 

0 . 

.. 0 \ 

( \ 

w = 

0 


" 1 , and A = 

A2 


\ 0 

0 . 

•• WN J 

\ Aat / 


Then the solution is expressed as follows, 

X = C’^A^w^v, (30) 

where C = A^w^A. And the value at boundary point {xb, yb) is evaluated by the value of the unconstrained 
polynomial P*{xb,yb)- 

Using this interpolation process, the intermediate boundary variable is calculated. Then, to determine 
the pressure on the boundary, a second intermediate boundary density and pressure are calculated by 


Pi* = P6(l + 7mMa2)^ 
Pi* = plil + T^Mal)^ 

P*b* = pI 
P*b* = pI 


for leeward surface where(U 5 — Vwaii) • n > 0, 


for the other situation, 


(31) 

(32) 


where n denotes the normal direction of a boundary surface, and Man denotes the Mach number defined by 
the ratio of the orthogonal relative velocity to the boundary and the local speed of sound. The density and 
pressure assigned at the leeward surface are extremely important for supersonic flow, since the extrapolation 
from the downstream to the leeward surface is physically invalid at supersonic speed. Therefore, the stagna¬ 
tion point condition for one dimensional isentropic flow is adopted to construct the upstream flow condition 
at the leeward surface. Obviously, the intermediate boundary velocity is not the no-slip boundary condition 
or no-penetration boundary condition. Then, W^* is transformed to Wf, to satisfy the boundary condition. 
To implement the isothermal no-slip boundary condition, we also assume that the pressure derived from 
Wb is identical to the pressure derived from W^*. And the velocity and the temperature are replaced by 
the boundary velocity and temperature respectively. That is. 


** 

Pb=Pb . 

(33) 

Ufo — ^walh 

(34) 

Pb=P*b*/iRT^au), 

(35) 

jp 1 tt2 , + 2 ** 

PbEb — ^PbVb + 2 ■ 

(36) 

The Euler slip boundary condition (reflection boundary condition) presents. 


Pb =P*b*y 

(37) 

\Jb = {n- U^a«)n + U,* - (n • U^n, 

(38) 

Pb = P*b*, 

(39) 

PbEh — ^Pb^h + 2 * 

(40) 

The no-slip adiabatic boundary condition gives, 


Vh =pT^ 

(41) 

wall: 

(42) 

** 

Pb = P , 

(43) 

1 2 K -\- 2 

PbEb — ^PbE^b + 2 ' 

(44) 



In this study, the solid boundary is stationary, thus \Jwaii = 0. Then, the boundary variables Wb at several 
boundary points {ixbj,ybj)\j = 1, < dim(X)} must satisfy the constraint. 

Constraint : Af,X = v;,, (45) 


where Vf, = {vb^i,Vb^ 2 , '^b,j denotes the every single component of 'Wbj, and 


Combining with Eq. 


Ah = 


/ 1 - 
1 - 


a^°,2%\2 


1 

•^ 6,1 yb,l 

^rn-1 1 

^b,2 yb,2 


<i2/m \ 


XbaVba 


\ 1 ... ^b,Myb,M 

the constrained solution is 
X = C-\I - Al C,-1A,C-1)A^w2v 
with Cb = AbC-^A[. 


VbM ^b^y^M ) 




b Vb, 


(46) 


(47) 

(48) 


So we fully determine the polynomial Pix, y) at the fluid side. If only with the consideration of one 
boundary point, then the current interpolation reproduces the least square method employed by Seo^. 
Furthermore, the unconstrained solution can be reused for calculating the constrained solution (for instance, 
C~^). Therefore, when the boundary variables cannot even be predetermined before the interpolation, the 
current interpolation method is identically efficient. In this study we only consider the linear polynomial 
interpolation. According to different boundary conditions, the derivatives at the interpolation points (or 
ghost points) on the solid side are determined by applying symmetric boundary condition. For the isothermal 
no-slip boundary. 


dp _ 9p. dp _ dp 

Qjjffiuid Isolid’ 1 fluid Isolid’ 

(49) 

9U„ _aU„ 9U„ _ aU„ 

1 fluid 1 solid’ 1 fluid Isolid’ 

(50) 

5Ut, aut, aut, _aut, 

1 fluid Isolid’ 1 fluid Isolid’ 

(51) 

d{pE) _d{pE) dipE) _d{pE) 

1 fluid 1 solid’ 1 fluid Isolid’ 

(52) 

for the Euler slip boundary, 


ap _ dp dp _ dp 

i^j^lfluid Isolid’ 1 fluid Isolid’ 

(53) 

au„ au„ au„ au„ 

Ifluid Isolid’ Ifluid g^ Isolid’ 

(54) 

aut aut aut au* 

g^ Ifluid g^ Isolid’ g^ Ifluid g^ Isolid’ 

(55) 

d{pE) _ d{pE) d{pE) _d{pE) 

g^ Ifluid g^ Isolid’ g^ Ifluid g^ Isolid’ 

(56) 

for the adiabatic no-slip boundary. 


dp. _ dp dp _ dp 

g^ Ifluid Isolid’ g^ Ifluid g^ Isolid’ 

(57) 

au„ au„ au„ au„ 

g^ Ifluid gj^ Isolid’ g^ Ifluid g^ Isolid’ 

(58) 

aut aut aut au* 

g^ Ifluid gj^ Isolid’ g^ Ifluid g^ Isolid’ 

(59) 

d{pE) _ d{pE) d{pE) _ d{pE) 

g^ Ifluid Isolid’ g^ Ifluid g^ Isolid 

(60) 


where ^ and ^ denote the partial derivatives along the normal direction and the tangential direction of 
boundary surface respectively. Up to now, the linear interpolation polynomial at the solid side is totally 
determined, so is the variable at the interpolation point. This interpolation procedure is adopted for both 
low speed and supersonic flow simulations. 
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FIG. 3. The mesh at left bottom corner (left figure). The time series of the total mass of all the fluid points for the 
lid-driven cavity flow at Re = 1000 (right figure). 


IV. NUMERICAL RESULTS 

This section includes numerical examples from the subsonic incompressible limit to the hypersonic com¬ 
pressible viscous flow computations. 


A. Subsonic 

The isothermal no-slip boundary condition is adopted for all the subsonic flow simulations. 


1. Lid-driven cavity flow 

The lid-driven cavity flow is a classical benchmark problem. With the Cartesian grid method, it is not easy 
to simulate internal flows, since the method theoretically cannot fully guarantee the mass conservation in the 
computational domain. Even though, it is still valuable to try the cavity case to estimate the conservative 
property of the Cartesian grid method. 

The computational region is defined on a domain of [—0.01,1.01] x [—0.01,1.01], where two vertical walls 
locate at x = 0 and x = 1, two horizontal walls locate at y = 0 and y = 1. The length of the cavity is 
Lq = 1. The upper wall is moving from left to right with a constant velocity. The isothermal boundary 
condition is adopted. The boundary values are given as follows, 

U=(Co,0), r=l, 0<x<l,y = l, (61) 

U = (0,0), T = 1, the other boundaries. (62) 

The number of grid points in each direction is 261. Since a margin adjacent to the solid wall is preserved, the 
internal area of cavity is approximately discretized by a 256 x 256 mesh (see the left bottom of the cavity in 
figure[3l). For the later simulations, the vertical and horizontal boundaries are cut by the grid lines similarly. 
The density is 1, and the temperature is 1 at the beginning. The figure |3| shows the total mass inside the 
cavity. The total mass should be a constant for a steady flow. However, the constant slope of the profile 
indicates that the total mass increases with time increasing. The gain or loss of the total mass depends on 
the relative position of the grid point to the boundary point and the specific boundary treatment. The mass 
changing rate is about 0.1% per unit time for this test case. The numerical results are compared with the 
Ghia’s reference value. The velocity profiles at time t = SOLq/Uq, when the total mass is not changing too 
much, are plotted in the figure S] 

Although the boundary condition is not fully conservative, the velocity profiles are still quit good. The 
viscous effect can be simulated correctly in this case. 




























10 




FIG. 4. The velocity U profile along the vertical central line and the velocity V profile along the horizontal central 
line for the lid-driven cavity flow at Re = 1000. The grid is 261 x 261 covering the range of [—0.01,1.01] x [—0.01, 1.01] 




FIG. 5. The mesh around the cylinder (left figure). The wake length vs. Reynolds number for the flow past a circular 
cylinder (right figure). 


2. Flow past a circular cylinder 


The two dimensional flow passing through a circular cylinder is considered to test Cartesian grid method 
for steady flows at very low Mach numbers. The free stream Mach is about 0.08, and Reynolds numbers 
based on the cylinder diameter, D, is Re = 10, 20, 30. A rectangular computational domain is used with the 
inlet and outlet vertical boundaries at a; = —IbD and x = 15Z1, and the two horizontal far field boundaries at 
y = ±1511, respectively. The origin coincides with the center of the cylinder. Computations are performed 
using a grid of 903 x 903 points. The grid size is about 0.033Z1. The curved boundaries are cut by the grid 
lines as shown in the figure [5] on the left. The figure on the right shows the wake length at the rear of the 
cylinder. The pressure coefficient is defined by Cp = {p — Po)/{1/2{pqUq)), where the subscript ’0’ denotes 
the upstream or inlet flow condition. The reference data is extracted from the references^^Ti^i. Figure [S] shows 
the pressure coefficient along the cylinder surface. The horizontal axis represents the angle of circumference 
where 9 = 0 and 9 = 180 correspond to the stagnation point and the rear of cylinder respectively. The 
agreement is quite satisfactory. It should be emphasized that the flow quantities at solid boundary are quit 
smooth in current computation owing to the constrained weighted least square interpolation. 
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FIG. 6. The pressure coefficient along the circular cylinder surface. 


grid size (da;) 

Li error 

Loo error 

0.03333 

0.01664 

0.00831 

6.781 X lO”*" 

2.346 X 10"® 
8.115 X 10"® 

2.140 X 10“* 

7.302 X 10"® 
4.697 X 10"® 


TABLE I. The pressure discrepancy between the upper and lower surfaces of the cylinder for the symmetric flow on 
the asymmetric distributed meshes at Re = 20. 


3. Boundary accuracy 


The low Reynolds number flow passing through a circular cylinder is planar symmetric in previous example. 
Therefore, the flow quantities should be exactly identical from the upper and bottom half cylinder surfaces. 
However, if the cylinder is not placed symmetrically (mirror reflection relative to the x-axis), such as a small 
shift of the cylinder vertically, the geometrical locations of the upper and lower half cylinder surfaces will 
not be symmetric. The error introduced from this kind of discrepancy can be considered as a measure of 
the accuracy of the Cartesian grid method. The Li error is defined as Erri = i \P{ff) — P{2 tt — 0)\d9\ 
the Loo error is defined as Ervoo = maxe{P(6*) — P{2'k — 6*)}. We calculate the same problem on three 
successive refined meshes, i.e., 301 x 301, 602 x 602, 1204 x 1204. The computational domain covers the 
area, [—5.004,4.996] x [—5.01,4.99], and the center of the cylinder is (0,0). Obviously, the origin and the 
center of cylinder are mismatched. The flow condition is the same as the setting for the Reynolds 20 case 
in the last example. The errors are shown on the table [H and also plotted in the figure 0 The current 
Cartesian grid method achieves an overall second order accuracy, and has only first order accuracy at some 
irregular points. This example totally excludes the influence of the flow solver. Thus, it is a good example 
for the accuracy evaluation for any Cartesian grid method. 


4. Flow over a square cylinder 

The 2D laminar flow around a square cylinder with length D mounted at the center of the channel is 
investigated. The configuration is the same as the reference^^. The blockage ratio is fixed at B = 1/8. So 
the height of channel is 8D. In order to reduce the influence of inflow and outflow boundary conditions, the 
length of the channel is set to be L/D = 35. The square cylinder locates at the length of lOD from the 
inflow boundary. The grid is 1416 x 341 points. The incoming flow’s velocity profile is given as a parabolic 
curve. The maximum velocity in the middle of channel is 0.1. The initial gas density is uniformly 1. The 
inlet pressure is 1, and the pressure of outlet is given based on a typical pressure drop for tube flow. Figure 
[8| shows the pressure coefficient contour at Re = 20 and the time series of the unsteady vertical velocity 
at location (3.0,0.0) when Re = 100. The Strouhal number derived from current method is about 0.137, 
which is identical to the reference datai^^. The wake length behind square is compared with reference data 
in figure [HI It also coincides with the numerical simulation by Breuer— . 
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FIG. 7. The pressure discrepancy between the upper and lower surfaces of the cylinder vs. grid size on log scale. 




FIG. 8. The left figure is the pressure coefficient contour when Re = 20. The right figure is the vertical velocity 
series at the location (3.0, 0.0) behind the square when Re = 100. The corresponding Strouhal number is 0.137. 



FIG. 9. The wake length vs. Reynolds number for the flow past a square cylinder. 
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FIG. 10. The contour of the pressure coefficient (left) and the pressure coefficient (right) along the cylinder surface 
for the flow around a circular cylinder at Mach = 3. 



X 


FIG. 11. The pressure coefficient contour and the stream lines for the flow past multiple objects when Ma = 3. 


B. Supersonic flows 

1. Supersonic flow around a circular cylinder 


Flow around a blunt body is a very important test case for supersonic flow solver. We use this case to 
check the robustness of the current Cartesian grid method. The gas initially has a velocity of Mach 3. The 
Reynolds number defined by upstream condition is 1420. And the upstream pressure is 1. The grid with 
363 X 303 mesh points covers the range of [—4.0, 8.0] x [—5.0, 5.0]. Isothermal no-slip boundary condition is 
employed on the cylinder surface = !)• The scheme survives initially at the strong rarefaction wave at 

the rear part of the cylinder. The figure [TUI shows the pressure coefficient contour and the pressure coefficient 
along the cylinder surface at the steady state. Then we put another two different geometries, say, a plate and 
a triangle, into the computational field. The flow condition is identical to the previous simulation. Figure 
m shows the pressure coefficient contour and the stream lines. This simulation demonstrates the robustness 
and flexibility of current Cartesian grid method. 
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FIG. 12. The left figure illustrates the mesh around the step. The right top figure shows the density contour for 
Mach 3 step flow. The right bottom figure is the contour of p/p^. The grid is 304 x 102 covering the range of 
[0.0,3] X [-0.05,1.05]. 

2. Mach 3 forward step flow 

This two dimensional test problem has been studied by Chiang et al^ using Cartesian grid method. 
However, in their figures, the grid lines are superposed on the boundary. In our simulation, the computational 
domain is [0.0, 3.0] x [—0.05,1.05]. Two horizontal boundaries are placed at y = 0.0 and y = 1.0 respectively. 
A step with height 0.2 is located at x = 0.6. The upstream velocity is {U,V) = (3,0). The gas in the 
computational domain is uniform at the beginning, whose density is 1.4, pressure is 1.0, and velocity is 3. 
The Euler slip boundary condition is imposed at all solid boundaries. The horizontal direction is discretized 
by 304 points. Then we change the number of the vertical discrete points. The length of Mach stem at the 
top surface varies when the vertical discrete point number changes from 100 to 110. This is because that 
the minimum distance to the top boundary changes dramatically when in the vertical discrete point number 
changes. And if double the mesh points in each direction, the numerical solution becomes insensitive to the 
small changes of the vertical discrete point number. Finally, the solution indeed converges to the reference 
solution^ when the mesh is fine enough even if the shock exists near the solid boundary. The density 
distribution in the wind tunnel and the entropy (p/p^) at time t = 4 are shown in figure [T^ No special 
treatment is adopted at the corner of the step. 


3. Double Mach reflection 

A wedge is placed at the origin with 30 degree angle to x direction. The computational domain is 
[—0.2, 3.464] X [0.0, 3.0]. The gas at x > 0, is at rest with density 1.4, and pressure 1. The flow condition at 
X < 0 is as follows, p = 8,U = 8.25,p = 116.5. A mesh of 900 x 730 points is used. The Euler slip boundary 
condition is imposed at the solid boundaries. The figure [13] displays the density contours. This case has 
already been studied by Cartesian grid method a^^d^ . In Forrer’s worbH, the whole structure moves slower 
than benchmark results. Our simulation results coincide with the results obtained by Pember— and the 
results from the calculation of a conformal mesb^. 


4. Viscous shock tube 

This is a viscous problem introduced by Daru and Tenau d^^d'^ . Sjogreen and Yee^ and many other 
researchers also studied this problem. The gas with different thermal condition is initially at rest on both 
sides of X = 0.5, and separated by a membrane. Then, the membrane is removed and wave interaction 
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FIG. 13. The mesh is 900 x 730 covering the range of [0.0, 3] x [—0.2, 3.464|. The bottom left figure shows a portion 
of the mesh. The bottom right figure shows the density contour of double Mach reflection. 



>- 



FIG. 14. The left hgure shows the mesh around the right bottom corner of the computational domain. The right 
figure shows the density contour of viscous shock tube. The grid is 2043 x 1023 covering the range of [—0.51, 0.51] x 
[-0.01,0.51] 


occurs. The compressible Navier-Stokes equations with adiabatic no-slip boundary conditions are imposed. 
The solution develops complex two dimension shock/shear/boundary-layer interactions, which depend on 
the Reynolds number. The dimensionless initial states is given as follows. 


Pi = 120 , Pi = 120/7, p, 


1.2,Pr = 1.2/7, 
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where subscripts ”1” and ”r” denotes the left side and right side of a; = 0.5 respectively. The Prandtl 
number is 0.73. The case of Reynolds number 1000 to time t = 1.0 is simulated. The computational domain 
is [—0.51,0.51] X [—0.01,0.51]. A box of [—0.5, 0.5] x [0.0,0.5] embraces the fluid domain. 2043 x 1023 
grid points are used in this simulation. The bottom and two lateral boundaries are adiabatic with non¬ 
slip condition. The top boundary is reflection boundary (Euler slip boundary). Generally speaking, the 
structures of this problem (figure I14|l are quit similar to that derived in reference o^^’^^ . This simulation 
demonstrates that current Cartesian grid method is capable of simulating complex supersonic viscous flow. 


V. CONCLUSION 

In this paper, we present a Cartesian grid method for complex immersed boundary problems, and pro¬ 
pose a simplified gas kinetic scheme for subsonic and supersonic flow simulation. A constrained weighted 
least square method is employed to update the physical quantities at the interpolation points. Different 
boundary conditions, including isothermal, adiabatic, and Euler slip conditions, are presented by different 
interpolation strategies. The new method is capable of simulating inviscid and viscous, compressible and 
near incompressible flow problems. The numerical results demonstrate that the new method converges to 
the correct physical solutions for both subsonic and supersonic flows as the mesh refines. The current scheme 
is robust under various flow condition from low Mach number to high Mach number flows, from inviscid 
to viscous ones. All these test cases are calculated by the same Cartesian grid method. The interpolation 
procedure proposed in this study provides a smooth distribution of physical quantities at solid boundary 
and can tackle with arbitrary geometry. The current methodology can be further extended to the other flow 
systems, such as to the numerical schemes^^— for both continuum and rarefied flow computations. 
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